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ABSIT RACI 


Three forms of the airplane spin equations of motion, derived by 
Buchler in Reference [1], form the basis for the development of a 
computer program designed to seek dynamically stable equilibrium 
solutions of a spinning aircraft. The program incorporates two 
solution techniques: one based upon Euler integration, the other, a 
version of minimization by gradient search. Secondary programs are 
developed to (1) generate ЕТЕР: glide parameters for use in the 
validation of the equations of motion, and (2) evaluate equation residuals 
obtained froma grid of initial conditions over the potential solution 
space. F-111l and F-4 aerodynamic force and moment models were 
utilized to evaluate the solution methods and equations of motion. The 
numerical results indicate that the F-111 and F-4 data are not repre- 
sentative of the actual aircraft and, therefore, it is highly unlikely 
that dynamically stable equilibrium solutions can be achieved from 
these models. The utility of the two solution methods is evaluated and 
the numerical results are analyzed in order to gain insight into the 
optimal application of the three forms of the equations of motion. The 
paper concludes with a discussion concerning the qualitative validation 


of the equations of motion. 
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TABLE OF SYMBOLS 


The definitions of the symbols used in this paper are as follows: 
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Un 


Center of gravity 


Mean aerodynamic chord, feet 


.Acceleration of gravity, ft/se ес 


Airplane principal moments of inertia, slug-ft 
Airplane mass, slugs 

Load factor; nugaber of "2's 

Spin Radius, feet 


Criterion function value; equal to the sum of the 
absolute values of the equation residuals, SUM 


Airplane wing area, Te 

lime, seconds 

Increment of time, seconds 

Airplane cg velocity, ft/sec 

Airplane principal axis system 

Reference Cartesian axis system 

Airplane altitude, feet 

Rate of ascent, ft/sec 

Airplane angle of attack, alpha 

Direction cosines in terms of Euler are (Table I) 
Airplane sideslip angle, Beta 

Spin rate, per second 

Aileron, rudder and elevator deflections, degrees 


Airplane principal axis system orientation Euler angles 
with respect to the X, Yi, Z1 system 





о Y, ф. Airplane orientation in pitch, roll and yaw with respect 
P,r, Y to the x Yi, Ay system 


in Air density, dn 


W Aircraft weight, lbs 
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Ia INTRODUCTION 


A. BACKGROUND 

Aircraft flight tests frequently reveal significant differences between 
the predicted and the actual spin characteristics of high performance 
systems. Such differences often result in expensive engineering changes 
to the aircraft designeand tend to reduce the mission capability of the 
aircraft. Furthermore, if the final spin characteristics are such that 
pilots are not permitted to practice the recovery from intentional 
spins, then many will consider the outer regions of the performance 
envelope with trepidation resulting in further degradation of the weapon 
system's effectiveness. The need therefore, to accurately predict the 
spin characteristics of proposed aircraft designs is of paramount 
importance. 

Historically, efforts towards spin performance prediction have been 
concentrated in two areas: the development of system time histories 
and the analysis of free-spinning model test results. System time 
histories involve the simultaneous solution of the aircraft equations of 
motion with the required force and moment coefficients being deter- 
mined from wind tunnel model data. Solutions are obtained using 
digital computer programs which numerically integrate the appropriate 
equations of motion. There has been minimal success towards 
achieving a time history method that will accurately predict spin 
performance due to the critical dependence upon entry conditions. 
Free-spinning model tests utilize instrumented dynamically scaled 
models which are generally deployed into vertical wind tunnels. This 


method however has a limited capability in simulating the full scale 
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spin environment (Reynolds, Mach, and Froude Numbers and the spin 
entry mechanism) and often the tests fail to identify all of the possible 


spin modes that may be experienced by the full scale vehicle. 


Be GENERAL 

This paper represents the final report on the development of a 
digital computer program designed to analytically predict the airplane 
«Spin performance using a unique set of aircraft equations of motion 
and simplifying assumptions regarding the nature of spin mode 
solutions. The primary hypothesis upon ED the development of 
the computer — was predicated was the unique dynamically 
stable equilibrium (steady-state) solution for each aircraft spin mode. 
This assumption implies that all motion will ultimately be damped to 
a steady-state solution. Using the above hypothesis, the dynamically 
stable equilibrium solution may be determined directly without the 
encumbrances ofa time history approach, with its required dependence 
upon entry conditions. The hypothesis ignores, however, the possible 
(and perhaps most probable) existence of a dynamically unstable 
equilibrium solution which could also yield a stable trajectory in 


state space', 


SS COPE 


A rigorous research effort into the development of a means of 


de 


predicting airplane spin modes would involve the following stages: 


J. Literature search 
[1. Derivation of the mathematical model 
lil. Verification of the mathematical model by using 


time history comparisons with baseline data 


l. Development of a computer program designed to 
seek spin solutions 


12 





V. Verification of the computer program by utilizing 
wind tunnel data modified for simulator use 


VI. Determine and/or predict "real world" spin modes. 
Stages Iand II are described in References [7] and [1]; while the 
base line data required for phase III is included in Reference [2]. The 
purpose of this effort was an attempt to accomplish steps IV and VI, 
This rather large undertaking was conducted concurrently with LCDR 
‘Buehler, who developed the applicable equations and coordinate systems 
required by stage Il. The limited time available precluded a numerical 


validation of the computer program (phase III). 


13 





il, ANALYTIC CONSIDERS TICHA 


Poe OORDINATE SYSTEM 

The coordinate system used in the mathematical model is a com- 
bination of the airplane principal axis system and a fixed cylindrical 
system. The vertical (z) axis of the fixed cylindrical system repre- 
sents the initial central axis of motion. The cylindrical system 
locates the aircraft center of gravity in terms of the altitude coordinate 
(20) the spin radius coordinate (R) and the angular position coordinate 
(BY), The orientation of the aircraft with respect to the cylindrical 
system is described in terms of Euler angles. A cartesian coordinate 
system is fixed at the cg position on the R vector with its axis (ху) іп 
the (-R) direction and the (21) axis is oriented іп the (-z) direction. 
The cartesian system provides a reference for the orientation Euler 
angles whereby zero values of (8), (W), and (6) would yield an upright, 
wings level airplane with its principal (x) axis pointed inward in the 
(-R) direction. Figure l depicts the details of the coordinate system; 
the positive directions of the position coordinates are also shown. 
As Buehler noted "The rationale behind this choice of coordinate 
system is that it more simply (in a mathematical sense) represents 
the motion being modeled.'' [Reference 1]. The advantage of this 
hybrid coordinate system is that it obviates a steady-state spin 
solution. The motion of the spinning aircraft is easily recognized by 
constant values of the six independent coordinate variables (R, 20, N, 
e, Y , Ó ) Additionally the complexities of the stall, departure, and 
incipient spin motion have no impact on the results since the objective 


is to determine the steady-state spin modes. 
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Fig. ). The Coordinate System Used to Model the Spin 
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Fig. 2. A General Schematic of Fuler Angle Relatie nmp: 
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EL EULER ORIENTATION ANG. 

The Euler Angle relationships utilized in the mathematical model 
involve rotations which are measured relative to the vertical (21) axis 
and to the intersection of the (x, y) plane and the (x4; y)) plane. This 


system yields equations for the direction cosines as listed in Table I. 


able TT: - Direction Cosines in Terms of Euler ongles 


Ug PC TE bas n (On AA pt ER. jT PD nm s e 


Cosines of Angles between X, Y,Z and = M 
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сц = cos $ COS y 
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~ sing siny 






а = cosgsiny к ү 
15 оз) = —singcosy 


+ sing cos y cos 8 + cos 9 cos cos 6 


сз sin é sin é саз sin é cos ó 834 — 00398 


Inherent in this orientation system, however, is the difficulty in 
visualizing the final orientation without ss o m sketches or digital 
graphics equipment. This can be readily appreciated by referring to 
the general schematic of the Euler angle relationships depicted in 
Eure 2. In order to alleviate the visualization problem, Buehler, 
in Reference [1], derived relationships between a set of orientation 
angles with a sequence of pitch M roll (W „) „and then yaw (9 x 
and the Euler angles depicted in Figure 2. A The equations for 
А reafter, and in the computer program, the orientation angles will 


be referred to as the 'ordered' Euler angles, while the orientation 


angles depicted in Figure 2 will be referred to as the ‘reference’ 
Euler angles. 
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converting the ordered Euler angles to the reference Euler angles are: 


| | Ə 
Ó -y- Tao | Too (Se) гл). ara) (1) 


Y = 7 ы Ti ar (e) Tau (%)} (2) 


PR 2) 


заето (Св) Ta w (£ y) 


Smf Nus ар ES 1] / An (E) 


Ө - 2 Tow” 


The equations for obtaining (е), у) and (9, from the reference 


Euler angles are: 


analy tee А По (гаа 
у= ф + 7м — 3 | efe (2)5 w(t = 2) (4) 


S 1 22 +6) Cos 12 +6 


V) ^0. I os 12 -Ө) _ Ze Y — 
Op = Ar а — Tam” mese, (5) 


Cos (22-9 2) Cos (— a = ©) 
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S/N suf rw | Zk ran! | Towle, Tan) sn FT ale zs / 
LET 
Q, = Z - 2⁄2" XI c LS 
5и 


uf roe fre ее cos ү == D 
"us 
Jf LST ше 


The above equations have been numerically verified and validated 
for quadrants І апа ІУ (+ 90): This is sufficient for an aircraft 
oriented in an upright spin, however, for the case of an inverted spin 
where the aircraft will encounter quadrants П апа ПІ, the above 


relationships are invalid. 


ЕНЕ БООАТІОМ5 ОЕ MOTION 

The spin is an uncontrolled large angle, six degree of freedom 
motion experienced by an aircraft in the stalled aerodynamic region. 
Since the spinning phemomenon is a complex motion which is influenced 
by а host of nonlinear variables, the equations which describe that 
motion are both nonlinear and coupled. Buehler, in Reference [1], 
derived the aircraft spin equations of motion utilized in the computer 
Er ram. The equations and modifications thereto are described 
below. 

l. The Full Equations - 

The equations of motion were cast in a more useable form by 
solving each equation for the corresponding acceleration parameter. 
This form avoids obscuration of the physical nature of the problem and 
places the numerical analyst in a better position to monitor, trouble- 


shoot and interpret iterative solution results. 
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2. The Modified Equations 
The full equations were further modified by recognizing that 
only dynamically stable (or quasi-stable) equilibrium solutions were 
sought. Thus angular orientation rates (8,% Р Фф ) меге equated to zero, 
which in turn eliminated a considerable number of terms. The R 
parameter was retained because it effects the values of &, le and DE 
to the extent that when ignored, it seriously impeded the iterative 
solution efíort. 
3. The Short Equations 
The modified equations were further simplified by equating the 
R parameter to zero. This form is useful when seeking purely steady 
state solutions using the more mathematically elegant solution methods 
such as the gradient search method described in Section H. 


The full, modified, and short equations are listed in Appendix I. 


D, ANCILLIARY RELATIONSHIPS AND DEPENDENT VARIABLES 
In order to utilize conventional aerodynamic force and moment 
coefficients in the equations of motion, certain ancilliary relationships 


were needed. As derived by Buehler, and as depicted in Figure 3: 


Angle of Attack: 


A = Jay | eu өзг OK T oom Zo (7) 
ou AX +A ОА * os X, 


Sideslip Angle: Р 


{= TA Le A 72 + 23 Z, (8) 
oly E +A sas AÑ Z e, = 


l. 





Relative Wind Velocity: 


2 AZ La Z (9) 
Ve7 1 R + (AR) x 
Roll Rate: 
os (4 + 7) ез +6 Cos 0 (10) 
Pitch Rate: 
ў. = (v +) oz - Sin È (11) 
Yaw Rate: 
д = (4 +) 4з + Ф (12) 
Á Y 
' 
' 
t 
v | 








Figure 3 


Definition of Positive& and P Angles 
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E AERODYNAMIC FORCE AND MOMENT СОЕ ЭСЕ 
The general form of the aerodynamic force and moment coefficients 


employed in the equations of motion are as follows: 


Cre = Cre, + Ce, Fle) + Ce РО) * Cz fco (13) 
5 Cr, кені) a (бə fup. * Се), “с бе) 


АЕ 


Cm, = Cm + Cm, EA + Cm, fg) + Cm, Fr) 
A CZ, TES) * Cu, f ($.) 2 En, / (je) 
xy 


(14) 


where a ie m represent the angular deflections of the ailerons, 
rubber and elevator respectively. 

All of the above coefficients on the right side of the equality sign 
are aerodynamic force and moment derivatives as defined in Table II. 
They must be experimentally determined and tabulated for a wide 
range of © and @ vlaues. f(p,q, or r) are functions of roll, pitch 
and yaw rates. 

The coefficients in Table II utilize € as a common characteristic 
length, thus allowing the use of a single length parameter input to the 
computer program. Since X and д are zero for the steady spin, 
aerodynamic derivatives with respect to these two parameters have not 


been included in Table II. 
СА 





SUBSCRIPT j 





Table II 


Definition of the Aerodynamic Derivatives 


The validity of any spin solution results is dependent upon how well 
the aerodynamic forces and moments have been modeled. One needs 
only to visualize the different dynamic environment that exists between 
the static model in a wind tunnel and the actual airplane expe riencing 
three dimensional stalled flow ina spin. Herein lies the greatest 
weakness in the analytic spin solution effort. The static wind tunnel 
data utilized to determine the above 42 coefficients provides, at best, 
a very crude model of the actual conditions encounte red by à spinning 
aircraft. However, for want' of better data, the conventional data was 
utilized in order to exercise the program and to seek a first-order 


approximation to à spin solution. 


Aerodynamic derivatives with respect toó and B have not been 
included sinceX% and are zero for the steady spin. 
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J. SOLVING THE EQUATIONS OFM IO. 

A casual examination of the equations of motion found in Appendix 
II reveals that they have the following characteristics: 

l. ordinary second order differential equations 

2. extensively coupled 

3, highly nonlinear 


4. force and moment terms are a function of three-dimensional 
flow over the geometry of the airplane. 


The above characteristics preclude direct analytical integration and 
thus one is forced to seek a numerical solution scheme which will pro- 
vide the greatest insight into the physical problem and also provide 
some degree of uniform convergence. Additionally, from numerical 
considerations the extreme length of the equations dictate the need for 
a method which maintains numerical significance and does not require 
excessive computational time. 

Two iterative solution methods were selected; the integration method 
which was utilized for steady and quasi-steady state PORTE and the 
gradient method for the purely steady-state solution efforts. In both 
instances, the set of equations could be visualized as a "black box'! 
function transform with independent variables as inputs and the equation 


residuals as outputs. 


fee trie INTEGRATION SOLUTION METHOD 

The integration method can be visualized as the incremental 
tracking of the trajectory of a multiple degree of freedom spring-mass- 
damper system. The incremental changes are generated by ''psuedo'' 
integration of the aerodynamic force nd moment inbalances manifested 


in the equation residuals. This set of incremental changes constitutes 
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the "trajectory". The psuedo integration is accomplished by taking 
advantage of the form of the equations whereby the residuals are 
essentially the accelerations of the six spin degrees of freedom. 
Application of a commonsncrement of time to elementary Ber 
integration formulae allows the computation of a dynamically consistent 
set of incremental changes in the independent variables and thus gener- 
ates new values for the subsequent iteration. Essentially, the method 
duplicates, in a numerical sense, what the aircraft experiences under 
actual flight conditions. 

| One of the major advantages of casting the equations of motion in 
a cylindrical reference frame, is the ease with which tumbling motion 
can be separated from the motion associated with a steady state spin. 


Constant values of Z R, $ ,0,V, and @ indicate a steady state 


0” 
solution, while non-zero values of Ө ) Y Mor indicate tumbling 
motion. The computational scheme discards any angular momentum 
inconsistent with a steady-state spin and thus on each iteration any 
Eubbng motion is suppressed, Such motion, if retained, could 
possibly carry the solution path (trajectory) through a steady state 
spin solution. 

The iterative procedure used in the computer was basically very 


Ro oe 


simple. In response toa set of initial conditions Gr 
Q ) and related constants, the corresponding accelerations were 
computed. Application of a common time increment (4t) to the basic 
Euler integration formulae yields a set of incremental changes which 
are then added to the previous set of independent variables to obtain a 


new set of initial conditions. The iteration process continues until the 


sum of the absolute values of the accelerations, denoted as S, is reduced 


24 





to less than some specified value. The '"spring-mass-damper' system 
thus moves closer to equilibrium as it converges to a steady state 
solution. 

The acceleration zn and R can be computed directly from the inde- 
pendent and dependent variables. The other accelerations ( Ж ; е Y 4 
Ф ) however, must be computed by a looping routine since the associ- 
ated equations are coupled to all of the angular accelerations. The 


looping ceases once a consistent set of accelerations is obtained. 


The incremental changes are obtained using the following equations: 


e e 


2 mou 


© (nr) Si) 


a ӘУ (15) 


Ф 


(Wn) 


Kae S A F(R any Art R&)h (17) 
М = N ‚(ш С (18) 


Owe) Ec (OA (19) 
Е-е, + Cu Me Oe a T 
Vind = Var +[YVat)e (21) 
Vir = Ya) ae T os де + 1р e£") e (22) 
Dam n Din + ( Ó At) € (28! 
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Figure 4 


A Schematic Representation of the 
Inte gration Solution Method 


26 





n «€ Ie 
Фо. a De £ ( | AS (b — )f (24) 


The constants (a,b,c,d,e, and f) are solution convergence factors 
which can be used to weigh the solution ''rate'' of a particular parameter. 
Also, one should note that the (©), (Y ), and (Ф) expressions are used 
only in conjunction with the full equations. A schematic representation 


of the iterative process is depicted in Figure 4. 


feet dE GRADIENT SOLUTION METHOD 

The existance of a dynamically stable equilibrium solution can be 
verified utilizing a set of equations which exclude all time dependent 
terms. A gradient solution method was developed to search for such 
solutions using the previously discussed short equations. The method 
utilized is commonly referred to as the optimal method of steepest 
descent. 

Conceptually, the solution space can be visualized as a series of 
peaks and valleys with solutions being the floor of any valley whose 
altitude is zero. The method is initiated at some 'geographic' position 
(initial conditions) where the "altitude'' is measured (criterion function 
value). If the 'altitude'' were non-zero, the slope (gradient) of the 
"hill" would be measured and then the independent variables would be 
incrementally ''marched' down the hill. The process continues until 


the independent variables reach the bottom of the ''valley"'. 


e [7], page 275. 
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Mathematically, the criterion function is defined as the sum of the 


absolute values of the equation residuals. A solution is obtained when 


7 
515150 (2.5) 


where 5: are the equation residuals computed from the given independent 
variable vector, X. (elo): 
The rapid reduction of S to zero is accomplished Dy iteratively 


choosing an increment $x; such that Š Š approaches -S where 





65 = > (I) Sx =, | (26) 


This is achieved by choosingóxas the product of the criterion gradient 


and a gain constant K; 


However, by noting that 


£ => (BR) (k FR) = -S (28) 


LEN 


N 


S 5326735 x (29) 


K can be expressed as 


NE re (30) 
Nis ШЕ 
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The subsequent value of the 5 element of the X vector becomes 


X, =X өз) 


{ AD 
f (A'r/) “(Q ) bx, (4) 


where the subscript n refers to a sequential iteration number. 


S is then computed as a function of X лігіне condition 


n+1 n+1 


EE (32) 


is satisfied, then the increment vector Sx, is reapplied, 


A X wr * dx, С 


and this 1-D search process is continued in the direction of the original 


gradient until 


К ОО И (34) 


whereupon new values of the gradients are computed. Since the process 
avoids the computation of the partial derivatives whenever the condition 
expressed in equation (32) is satisfied, significant amounts of computer 
computational time is saved and hence, the name 'optimal! gradient 
solution method is achieved. Figure 5 includes a flow chart of this 


technique. 
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MI, AIRCRAFT SPIN SOLUTION COMP UTER (nl 


NO GENERAL : 

This section presents briefly a general description of the digital 
computer program designed to solve the spin equations of motion for 
dynamically stable equilibrium solutions. A listing of the program as 
well as specific usage data and flow charts is presented in Appendix 
II. A program hierarchy is included as Figure 6 and an explanation of 


the program output is described by Figure 7. 


Pee ROGRAM USAGE 

The program reads the tabular force/moment coefficients and 
related aircraft constants and then commences to iteratively seek a 
solution based upon the first initial condition data set that is read. 
Upon termination of a particular solution, the program restarts by 


reading the next initial condition data set. 


@ZZINPUT DATA 
The user must provide the following information to the program: 

1. Type aircraft (F4, F111, etc.) 

2. Tabular values of the force and moment coefficients/ 
stability derivatives (as described in Table II) and the 
corresponding values of angle of attack and sideslip 
angle 

3. Aircraft mass 

4. Atmospheric density 

5. Mean aerodynamic chord 


6. Gravitational constant 


7. Total wing area 
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Airplane Spin Solution Program Hierarchy 


DE 





8. Principal mom€nts of inertia 
9. Control deflections 
10. Initial conditions of the independent variables 
The user also must stipulate the solution method, equation form, 
solution criteria, time increment, output format and solution con- 
vergence factors. 
The solution convergence factors available for the integration 


ате" 


method, allow the user to accelerate or decelerate the solution 
of any of the six independent variables. For example, by equating a 


particular factor to zero, that variable will remain constant, thus 


reducing the number of degrees of system freedom. 


PAR OUTPUT DATA 
The printed output format is similar for the integration and gradient 
methods. The following is printed: 
l. Listing of all input data under appropriate headings 
2. Iterative results (integration method) 
a. iteration number (М1) 


b. the number of looping iterations required to obtain 
a consistent set of angular accelerations (M2) 


c. the independent variables with their associated first 
and second derivatives 


d. relative wind velocity, angle of attack, and sideslip 
angle 


e. sum of the absolute values of the accelerations 
3. Iterative results (gradient method) 
a. iteration number (M1) 


b. the iteration number corresponding to the last 
gradient computation (M2) 
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c. the six independent variables with their corresponding 
residuals, incremental values, and partial derivatives 


d. relative wind velocity, angle of attack, and sideslip 
angle 


e. criterion value; the sum of the absolute values of the 
residuals 


f. K; the solution gain factor 
4. Summary output 
Upon norma! termination of the integration method, the iterative 
output is summarized. The listing includes the following data corres- 
ponding to every twenty-five iterations: 
l. iteration number 
2. independent variables 
3. dependent variables ((,/2 , Va ) 
Also, the above data corresponding to the smallest sum of the 
absolute values of the accelerations is printed. 
In the event of an abnormal termination, the program will print an 
appropriate message which informs the user as to the reason for the 


termination as well as the values of the offending variable(s). 


БКОСКАМ MODULE DESCRIPTION 

The program consists of nine modules; one main and eight sub- 
routines as described below: 

1. ` MAIN 

This routine is basically an input/output device. Tabulated 

data, fixed constants, and initial conditions are read into common 
memory, angular units are converted from degrees to radians; the 
optional output format is established and the first page is printed. 


The optional equation form and solution method, being determined, 
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a call is then made on the appropriate solution subroutine for sub- 
sequent computation. 
2. SOL NV Tm 

The integration solution method described in Section II is pros 
s irnred into this subroutine. The iterative results are printed and 
appropriate summary data is presented on the last page of the printout. 

5 GAWAIN 

This subroutine incorporates the optimized -gredient method-as 
EE ied in Section II. Its operation is similar to SOLVIT, however, 
no summary data is printed. 

4. EQNS 

Subroutine EQNS includes all of the independent and dependent 
equations. Calls are made primarily from either GAWAIN or SOLVIT. 
The current set of initial conditions is then used to compute the 
direction cosines, angle of attack and sideslip angle. Subroutine ANGLE 
Brea lcd to determine whether Alpha and Beta are within the appropri- 
ate definition limits. If so, the force and moment coefficients are then 
determined by interpolating the appropriate matrices using subroutines 
ieee R, SPLINI, and SPLICO. 

The residuals are computed using either the short, modified or 
long form of the equations. A looping routine is utilized with the latter 
forms to obtain a consistent set of angular accelerations. Finally, 
upon computation of the sum of the absolute values of the residuals, 
execution is returned to the calling program. 

О КЕПШЕ 

The equations of motion require that angular units be in radians 

and that the reference Euler orientation angles be utilized. Subroutine 


EULER is used to convert computer variables to the output form of 
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angular units of degrees and ordered orientation angles. The sub- 
routine is utilized by MAIN, GAWAIN and SOLVIT in order to allow 
output of the more usable ordered angles. 
ees NCE 
Subroutine ANGLE is utilized primarily by EQNS to check 
whether alpha or beta has exceeded the upper or lower limits of data 
definition. If the limits were exceeded, the iterative solution effort 
would be terminated and control returned to MAIN, ANGLE also reduces 
the absolute magnitude of angular displacements to less than 21i radians. 
foe Lx 
Subroutine FIX is used in conjunction with IBM supplied ERRSET 
to facilitate termination of an iterative solution in the event of an under- 
flow, overflow or divide check while control is in EQNS, Control is 
returned to MAIN rather than effecting program termination as is the 
usual case. 
INTER 
Subroutine INTER is utilized by EQNS to interpolate the tabular 
force and moment coefficient data for specified values of alpha and 
beta. INTER calls SPLINI for the cubic interpolation of each matrix 
СОМЕ Стог. 
DSESPLINI and SPLICO 
Subroutines SPLIN1 and SPLICO were written by Mr. Rod Kure 
and are maintained by the NPS Computer Center Staff. Their primary 


function is to provide cubic interpolation of a given matrix or vector. 
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V. AIRPLANE SPIN EQUATIONS PI ИО 


The airplane spin equations of motion were validated by establish- 
ing initial conditions for a power-off glide configuration and then 
applying them to the equations using the spin solution program. 

A small computer program, TRIM, was prepared to enable the 
selection of a set of dynamically stable initial conditions. The TRIM 
program flowchart, the program listing and the sample output is 


included in Appendix D. 


The condition of a power-off glide is mathematically described by 


three expressions: 


Tan (1-0) = А (34) 

W = хр о 005 (х-ө) + C, SiN (a -e&) P" 
— ( | 

e E mt (36) 


de 


Еге ч <,6,5.,. Ci, C and m are as depicted in Figure 8. 
The values for C and Ca are determined by transposing from the 


principal to the stability axis systems as indicated below: 


las а $e) Cos o. + (Cx "C ee Se) Sov ex (37) 


Co = (Cx * Cu. E +C + Ca, Se) Swa (38) 


Given constant values fora, w, P and S , one can solve equations 


(34)--(36) for the quantities $e , O , and V. Subsequently the variables 
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e 


Z, = V sw (a-8) (39) 
K = Vcos (2208) (40) 


can then be determined. 

The variables Lo R and e are thus established independent 
of the spin equations of motion but based upon the same aerodynamic 
force and moment model utilized by the spin solution program.  Appli- 
cation of these initial conditions to the spin solution program containing 
Buehler's spin equations resulted in an immediate convergence whereby 
the sum of the absolute values of accelerations (equation residuals) was 
equal to 0.0182. Table III provides a listing of the pertinent data 


utilized to achieve the equation validation. 


Table III 


Pertinent Data Relating to Equation Validation 


Д8 К Өр см, Т Se 
TRIM Input: F-111 force and moment model; air density = .1553 x ТІРЕ 
TRIM QUOI. IS 29 25855 -5 54. 
Output : = 
Integration 76. 99 ӘЗІР С 12.0 Secor 
Method 
Solution 


Notes: 1. Spin solution program imputas per Table VIII 
2. SUM = . 0188 
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The spin solution program requires special care when exercising 
it for straight and level flight. Due to EULER restrictions about 
accepting quadrant II and III and cardinal orientation angles, the 


initial condition values should be applied as in Table IV below. 


ACE TS I I 


Initial Conditions for Equation Validation 


En Qamqam s 00 Carr wam oo eae + s. AM Te ne í - accep ee ee ee 


VARIABLE VAIE 
P 5) 1% As used in TRIM Program 
Se Zo E, R As obtained from PDA Output 
Sa, Sr, ж 0.0 
4,0 ‚0001 
NSFLAG 0 


The results of the equation validation computer tests conclusively 
validate the spin equations of motion for the power-off glide configu- 
ration. This procedure can only be considered, at best, qualitative, 
for it does not test the equations ina time dependent environment and 
configuration. Complete validation would require the generation of 
time histories Bia subsequent comparison to base line data such as is 


included in Reference [2]. 
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V. NUMERICAL ANALYSIS AND CONCLUSIONS 


Pe ABERODYNAMIC FORCE AND MOMENT MODELS 

The primary purpose of this thesis effort was to develop an aircraft 
spin solution computer program. During the program development, there 
were no idealized data available with known spin solutions, therefore the 
author utilized actual static wind tunnel data to provide an aerodynamic 
force and moment model. This was done in one case with the objective 
of ultimately comparing actual aircraft spin performance with computer 
gene rated Есен уйне, 

Two models were utilized. Tabular F-111 data was obtained from 
Reference [3] and F-4b data were extracted from curves presented in 
Reference [4]. In both cases appropriate values were multiplied by 
b/c in order to normalize to the common dimension utilized in the 
computer program, i.e., the mean aerodynamic chord. The above 


data are included in Appendix E. 


SOLUTION METHOD EVALUATION 
1. Integration Solution Method 

The integration solution method was the primary vehicle 
utilized to develop the spin equations of motion. As discussed in 
Section II, the primary advantage lies with the characteristic of the 
equation residuals being equivalent to the appropriate acceleration 
terms. 

The method must be initiated with a set of initial conditions. 
Many sets of various initial conditions were utilized during the 
development but for the sake of brevity only the results from one 


arbitrary set will be discussed. 
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It was evident early in the development phase that the modified 
equations provided the greatest potential. The full equations were 
suitable only for verifying the stability of solutions previously obtained 
from the modified equations. When the full equations were used to 
evaluate a set of arbitrary initial conditions, the angular orientation 
rates rapidly increased to disproportionate values and induced 
tumbling motion which, in turn, resulted in a divergent solution. The 
modified equations, however, enable quasi-steady solutions and are 
responsive to small changes in the time increment and control 
deflections. 

The solution paths generated by the integration method were 
characterized by an initially rapid convergence to a minimum value of 
the criterion function. The orientation values reorient to a nearly 
optimal configuration; the sink rate seeks the acceleration minimum; 
and the spin rate and radius effect trade-offs to achieve a consistent 
velocity value. Occasionally, the solution path became divergent in 
that the program was unable to find a consistent set of coupled 
accelerations. This was evidenced by the printout of the number of 
iterations required by the acceleration ''do loop" within Subroutine 
EQNS. The 'do loop" limit is set at twenty iterations t o minimize 
computer time, however, only two iterations are usually needed for 
a stable solution path. Divergence is characterized by large and 
erratic changes in the independent variables such that alpha or beta 
exceed their respective fields of definition, thus resulting in solution 
termination. Solution path divergence can be minimized through use 
of a small time increment. The author determined that a range 
between .001 and .05 seconds proved optimal. The solution conver- 


gence factors were not utilized sufficiently to merit discussion. 
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Once the initial convergence was achieved, the individual 
solution paths seemed to vary ina seemingly arbitrary manner. 

While one variable tended toward a minimum acceleration value, other 
variables were producing increasing values. The variables continued 
to "wander'; changing magnitude and sign until either computer time 
was exhausted or the alpha or beta limits were exceeded. 

All solution efforts using F-4 data with pro-spin control 
“deflections resulted in solution paths which excecded tbe sideslip.angle 
limits. Inclusion of rotary balance data to compensate for differential 
sideslip angles along the fuselage might inhibit such yawing action. 

Utilization of the integration method requires that the user 
learn to iterate the output to optimize the achievement of a final 
solution. This process is as described below: 

a. Anarbitrary set of initial conditions is applied to the program 
to initiate a 500 iteration solution effort. This set of initial conditions 
should be such so as to approximate the anticipated spin solution mode. 

b. The initial condition set which yields the best convergence 
is modified by substituting the value for 2” corresponding to the 
minimum Z in the 500 iteration output. An alternative is to substitute 
the set of three orientation angles corresponding to the minimum 
absolute sum of the orientation angle accelerations. Similarly, R and 

$ can be substituted. The resultant set of revised initial conditions 
is then used to initiate subsequent "е efforts of shorter (150 
iterations) duration. 

c. The process is continued until a best solution is obtained. 

d. The final solution is tested for stability by using the full 


equations with a small time increment. 
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Figure 9 and Table V illustrate a typical solution effort using an 
arbitrary set of initial conditions. 
The integration method results lead to a nurnber of conclusions: 

a. The F-111 aerodynamic force and moment model pre- 
cluded finding stable equilibrium solutions. This can be attributed to 
any of four reasons: (1) there were no truly steady state solutions; 
only oscillatory modes exist; (2) the "right" set of initial conditions 
was never employed and therefore a steady state solution path was not 
intercepted; (3) the steady state modes were lightly damped and there- 
fore, may be very slow in convergence; the program wasn't exercised 
long enough to find the solution; and (4) the integration solution method 
itself precludes uniform convergence of all degrees of freedom. This 
would be due to self-excitation of the non-linear coupled acceleration 
te nmns. 

b. The F-4 data do not model the actual aircraft since all 
initial conditions yielded solution paths with divergent yaw rates. The 
actual aircraft does not exhibit this trait and, in fact, demonstrates 
several steady and quasi-steady spin modes. 

c. The method is cumbersome to use because it requires 
that the user iteratively adjust inputed initial conditions in order to 
achieve a solution. 

d. The radius rate terms incorporated in the modified 
equations enables the determination of quasi-steady state solutions 
using the F-111 data. 

a The orientation rate terms in the full equations were 


not compatible with the course integration techniques of this method. 
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Figure 9 


Schematic of a Typical Integration Solution Effort 
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f. The computer program needs to be modified to incorporate 
a subroutine which will compute the spin radius as well as the cylindrical 
inertial radius as in the present version. 
2. Gradient Solution Method 

The gradient solution method and its derivatives from optimal 
control theory have, perhaps, the greatest potential in the area of 
steady-state spin prediction. The method programmed by the autnor 
can only be considered crude when compared to:some of the more 
sophisticated methods which utilize hybrid applications. 

The results obtained from a typical F-111 gradient solution 
effort provide the basis for evaluation of this method as well as yield 
insight into the F-111 solution space. Table VI includes a summary 
of these results. One should note that there was an initial convergence 
of the angular variables in the first fifty iterations. After iteration 
#75, the solution stabilized but the criterion value alternated between 
27.0 and 102. 0; this was a result of the alternating residuals of the 
pitch and yaw equations. The sensitivity of the solution is extreme 
when one considers that a 0.4% change in the absolute sum of the 
independent variables results ina 2, 740% reduction in the absolute 
sum of the acceleration. Specifically, the alternating changes in the 
Gand @residuals are the result of incremental changes in the © 
and (D values of only home radians! 

The solution convergence becomes static due to the character- 
istically large gradients filo are and resultant small incremental 
step sizes. All gradient solutions attempted, resulted in no significant 
change in the values for 2. and R parameters. 


Two conclusions can be made from the gradient solution method 


results: 
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a. The method was unable to find a steady state solution; 
however, in one case, it was able to effect a reduction of the criterion 
function value by four orders of magnitude. 

b. The F-111 solution path is highly sensitive to small 
changes in the independent variable values. 
oes NALYSIS OF THE F-111 AND F-4 STEADY-STATE SOLUTION 

SEACHES 

The integration and gradient solution methods were both used to 
analyze the F-111 and F-4 steady state solution spaces. Initial 
conditions Were generated by a 'Criterion Function Search Program’. 
This program evaluated criterion function values from a course grid 
EN 350.624 sets of initial conditions for each aircraft. 

One can visualize the solution space as a geographical relief with 
the 'solutions! occurring at the valley floors of zero elevation. The 
search program essentially measured the 'altitudes' (criterion function 
values) in this area and found only a very few points whose elevations 
were less than 10*. Specifically, in the case of the F-4, only eleven 
sets of initial conditions generated criterion function values less than 
104 and only one set of F-111 initial conditions met this criteria. The 
numerical results of the criterion function search program are included 
in Table VII. 

These generated initial conditions were evaluated using the gradient 
and integration solution methods. The results are extremely signifi- 
cant because they reveal that it is highly unlikely that a dynamically 
stable equilibrium solution exists. These solution results are 
summarized in Figure 10 and Table VIII. 


- 


The two conclusions which can be made are: 
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1. The static wind tunnel data utilized to form the basis for the 
F-111 and F-4 solution spaces are of sucha nature as to preclude 
finding a dynamically stable equilibrium solution. There was no 
solution convergence toa criterion value of less than 25.0. 

2. The aerodynamic data utilized do not adequately describe the 
actual aerodynamic force and moment field experienced by the actual 
aircraft. This is evidenced by the failure to find even a quasi-steady 
solution similar to flight test results. 

The above implications are not surprising considering that the 
aerodynamics of the spin problem involve low speed, three-dimensional 
stalled flowand, as such, forces and moments are nonlinear. Thus 
the superposition of the separate effects of triaxial rotation, control 
deflection, and other static wind tunnel data cannot be completely 
justified analytically. The only wind tunnel test technique that 
presently appears to offer any possibility of yielding satisfactory 
aerodynamic data for the spinis the rotary balance force and moment 


measuring apparatus which is described in Reference [6]. 
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VI. SUMMARY OF THE MAJOR CONCLUSIONS 


The criterion function search revealed that, for the F-lll and F-4 
aerodynamic force and moment models utilized in this study, it was 
highly improbable that a dynamically stable equilibrium solution 
existed. 

The airplane spin equations described in Reference [1] were 
qualitatively validated. Quantitative verification is needed and it is 
recommended that the equations be cast ina form compatible witha 
time history analysis in order that the resultant data may be compared 
with the base-line spin data included in Reference [2]. 

The integration solution method is too cumbersome to be effectively 
utilized although it can provide information relating to quasi-steady 
solutions once the user has learned to iteratively manipulate the initial 
conditions. 

The gradient solution method demonstrated its potential by effecting 
Significant reductions in criterion function values to quasi-steady state 


solutions. 
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APPENDD A 


FULL SPIN BOUATIGNS 
The following equations constitute the full equations of motion for 


a spinning airplane as written in a cylindrical coordinate system. 


Z Equation 
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MODIFIED EQUATIONS 


Z o Equation: 
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SHOR MEOUATIONS 


2 о Equation: 
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APPEND 365 


AIRCRAFT SPIN SOLUTION FP ROGRAM 
A. GENERAL 
This appendix is intended for the Fortran programmer who is a 
potential user of the program. The included source listing contains 
.an.extensive introductory comment section with which the reader 
should familiarize himself prior to reading the other sections of this 


appendix. 


Eee tA INPUT 

The composition of the data deck is described in the source listing 
and also depicted in Figure D-1. One should particularly note the 
composition of data sets #2--#43. "The first card of each of the afore- 
mentioned sets ''tells'' the program which tabular data is to follow, 
its literal name, and whether the data set is a matrix, a vector, or a 
null set. The subsequent cards of each data set contain the tabular data 
except where no data is available (null set); in which case, the data set 
consists only of the first card. Each force/moment coefficient matrix 
is characterized by a unique integer assignment (NA) as specified in 
Table B-I. NA is the first value of each of the 42 coefficient data sets. 
There must be exactly 42 ''NA' cards however there is no restriction 
as to their sequence. 

Bach solution effort is initiated based upon the initial conditions 
and associated program variables provided by data set #45. An option 
is provided (integration method only) for the program variables to 


remain constant for subsequent solution efforts; only the initial 
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TABLE Dz I 


Force and Moment Coefficients, Stability Derivatives, 
and Associated 'NA' Assignments 
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conditions included as data set #48 need be provided. This option is 
effected by assigning NSFLAG = 2 or 3 as explained in the program 
comments listing and Table B-II. 

It is suggested that the user who wishes to test a large number of 
initial conditions assign LIMIT = 10, NSFLAG = 2 and then test all of 
the initial conditions on one job. This scheme will quickly reveal to the 
user which initial conditions exceed the alpha/beta limits without 
wasting valuable computation time. Subsequently, those initial con- 
ditions which appear to have a solution potential can be submitted with 
appropriate data set #45 values. 

Since every iteration generates a new set of initial conditions, the 
program can be restarted merely by including as data set #45 the appropri- 
ate values taken from the printed output. This option allows the user to 
check the stability of a quasi-steady spin solution by changing the time 
increment (DELT), one of the solution convergence factors (SCF), the 
control deflections (AILDEF, RUDDEF, ELEDEF), or the atmospheric 


density (RHO). 


RDA TA OUTPUT 

All output is formated for 130 character computer paper. An option 
is provided by NPRINT to suppress the iterative results; otherwise, all 
input data is echo checked, all significant iterative computations are 


printed and a summary of the integration solution path is provided. 


West ORAGE REQUIREMENTS 
The amount of storage required is primarily a function of the size 
of the CFORCM three dimensional matrix. The third index of the 


matrix is associated with, and is no larger than, the largest NA 


65 





TABLE B-II 


NSFLAG Integer Assignment 


NSE LAG SOLULION 
METHOD 
0 Integration 
1 Jntegration 
Z Integration 
3 Integration 
4 Gradient 
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EQUATION 
FORM 


Modified 
Full 
Modified 
Full 


Short 


SUBSEQUENT 
DATA SET 
COMPOSITION 


#45 
#45 
#48 
#48 


#45 





integer (NSMAX) corresponding to vector or matrix coefficient data sets. 
The first and second indices are associated with IMAX and JMAX, 

Prior to using the program the user should dimension A, B, and 
CFORCM in the COMMON/WORKA/ specification statements to 
A(IMAX), B(JMAX), and CFORCM(IMAX, JMAX, NAMAX). The 
WORKA statement must be identical in the following routines: MAIN, 
GAWAIN, SOLVIT, EQNS, INTER, and ANGLE, Common regions are 
utilized for most communications between subroutines and therefore 
care ms be exercised to ensure that the common statements are 


identical. 


EB, PROGRAM FLOW CHARTS 

Main and subroutine flow charts are included as Figures B-2 to B-9 
and are intended to be self-explanatory. The program source listing 
includes a sufficient number of comments to amplify the flowchart 
information. Subroutines SPLIN] and SPLICO were provided by the 
NPS Computer Center. No effort has been made to provide reference 
flowcharts as both routines have been extensively tested and validated 


by the NPS staff. 


F. PROGRAMMING CONSIDERATIONS 

One should note that calls on Euler are made in pairs. The angular 
variables are stored in the common memory, thus if the first conver- 
sion is from equation compatible to output nes the second call 
must effect the opposite conversion. 

Subroutine GAWAIN requires the computation of the partial derivative 


ofthe criterion with respect to a change in each independent variable. 


This is accomplished by incrementing a single independent variable by 
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FAC, The criterion is then computed using this value. The change in 
the criterion function value is then divided by FAC to obtain the partial 


derivative. 
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Figure B-6 


EULER Subroutine Flowchart 
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Figure B-7 


ANGLE Subroutine Flowchart 
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APPENDIX Ç 


CRITERION FUNCTION SEARCH PROGRAM 

The aircraft spin solution program was modified in order to evaluate 
the ''reasonableness'' of various sets of initial conditions. The result- 
ing Criterion Function Search program is designed to compute the 
criterion value for a given set of initial conditions. The initial con- 
ditions are generated by six nested 'do loops'' within MAIN. The 
dependent variables £,  , V and n are subsequently computed and 
compared to pre-established constraints; if the constraints are satisfied, 
the criterion is subsequently computed. If the value is less than IL the 
independent variables are printed. An abbreviated flowchart of the pro- 
gram is included as Figure C-l. Pertinent data relating to the initial 


condition 'do loops" 


and dependent variable constraints is included as 
Table C-I. The program is basically the same as the spin solution 
program except that SOLVIT and GAWAIN are not included and that 
MAIN and EQNS are modified to include the initial condition ''do loops'' 


and constraints. The range of the do loops is such that the initial 


conditions are all within a reasonable approximation to an erect spin. 
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[READ AERODYNAMIC MODEL DATA | 


SIX NESTED DO LOOPS > 





\ 

EX TO VARY INITIAL CONDITIONS / 
[COMPUTE DEPENDENT VARIABLES, 
CT DEPENDENT 

e ° VARIABLES : 


CONSTRAINT 


V LT. 


COMPUTE CRITERION FUNCTION VALUE 
_ (SUM) | 


GT. n 
== SUM 0 


ji. 





WRITE CURRENT INITIAL CONDITIONS | 


Figure C-1 


Abbreviated Flowchart for the Criterion Function 
Search Computer Program 





TABERE СЕТ 


Pertinent Data Regarding the Criterion Function 
Search Program Independent and Dependent Variables 


ima Condition 'Do Loop" Ranges 


Ez. R ж Gp Ur Ф 
-150. 0 25,0 10. 0 т TO - 5. 
-200. 0 100. 0 2070 -21.0 15. 0 -15. 
-250. 0 [50 30. 0 -31. 0 25.0 -25 
-300. 0 ¿5 0. 0 40. 0 -41. 0 35.0 -35. 
-350.0 325.0 50. 0 -51. 0 45.0 -45, 
-400. 0 400. 0 60.0 -61.0 D370 -55. 

-65. 
- 19, 
-85. 


QO OCO CO O O OPO < 


Depe adent Variable Constraints 


В = 0° + 40° 
&X = .1,0? to + 89, 0° 


Ү $ 589.0 (corresponds to 300 knots equivalent airspeed) 


IA 


Е 
n 2.0g's 


Fixed Constants 





Ё = . 001756 slugs/ft^ 
ЕТТІ F-4 
& = 100 21.0 
<A - 30. 0 
Pee -25.0 - 30.0 
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APPENDED 


TRIM PROGRAM 

















READ ¢, W, S, 
READ FORCE AND MOMENT VECTORS 
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SOMPUTE C. AMD 
EQNS (37) AND(38) | 
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COMPUTE 6 AND V 
EQNS (34) AND (35) | 
' COMPUTE Ż_ AND A 
EQNS (39) AND (40) 
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CONTINUE 


STOP) 


Subroutines required: SPLIN1 and SPLICO 
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APRENDO E 


AERODYNAMIC FORCE AND MOMENT MODETLS 


1. I 111 Data 
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slug a 


B. Attached data F-111 reduced to 6 x 9 


Wing Area 


Ti, F-4 Data 
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B. Attached data (F-4) reduced to 6 x 9: 
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